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Study guide 


You should schedule four study sessions for this unit, the main purpose of which is 
to help you revise for the exam. You will not need to use a computer. 


Sections 1 to 4 are of similar length and include revision activities. Section 5 is 
short and does not contain any activities. One possible study pattern is as follows. 


Study session 1: Section 1. 
Study session 2: Section 2. 
Study session 3: Section 3. 
Study session 4: Sections 4 and 5. 


Introduction 


You may have wondered what led us to choose the cover illustration for the course 
materials. In what sense, you might ask, does it illustrate ‘practical modern 
statistics’? Well, we must admit that a major factor in our choice was that the 
picture was free of charge ... But also, we felt it exemplified how unexpected 
patterns can be revealed by changing one’s viewpoint. In this case, literally: the 
picture is a satellite photograph of an alluvial fan between the Kunlun and Altun 
mountains in China’s Xinjiang province. 


At first sight, it is not at all clear what the photograph represents. Its strangely 
scale-free and organic look is reminiscent, perhaps, of a sample of biological 
material seen through a microscope, or an underwater picture of a jellyfish, or a 
strange nebula, or even a work of abstract art. What is very clear, however, is 
that this photograph of a portion of the Earth’s surface reveals an intricate 
pattern suggestive of a regular underlying structure. Most likely, this structure 
would not be immediately apparent from ground level. Modern statistics at its 
best is like that: it provides the viewpoints from which the structures and 
relationships in data that are otherwise concealed may be revealed. 


In this course, four topics of modern statistics have been introduced, each of 
which takes a different perspective, so as to bring out patterns in data. The 
techniques of medical statistics are used to reveal relationships between exposures 
and outcomes. Time series analysis is employed to tease out the underlying 
structures of data collected over time. Multivariate analysis can help make sense 
of large multi-dimensional data sets. And Bayesian statistics enables us to ground 
our inferences in expert knowledge. 


Another reason for choosing this particular cover illustration was that it 
represents a picture of part of our own planet. Statistics is about quantifying 
uncertainty and revealing patterns, not in any abstract sense, but as they relate 
to scientific questions about our environment and ourselves (hence the ‘practical’ 
in the title of the course). For this reason, throughout the course, real data and 
real problems have been used to motivate the methods described. 


In this unit, some of the methods discussed in the four books are used to address 
some practical issues related to a common theme of great importance: global 
climate change. The aim of the unit is to help you to begin revising for the 
examination. (The exercises at the end of each book may also be used for 
revision.) Not all the material required for the examination is revised here, and 
not all that is revised here will necessarily come up in the examination. 








The unit contains five sections. Sections 1 to 4 contain activities using ideas and 
techniques from Books 1 to 4, respectively. Section 5 is a brief description of some 
topics in medical statistics, time series, multivariate analysis and Bayesian 
statistics that have not been covered in M249. 


For more information, see 
NASA’s Earth Observatory 
website 

http: //earthobservatory.nasa.gov. 


Since the aim of this unit is to 
help with revision, and you will 
not have access to a computer in 
the examination, there is no 
computer work associated with 
this unit. However, you will 
need the Handbook and a 
calculator. 


1 The health consequences of climate change 


It has been estimated that, unless steps are taken to control the factors affecting 





climate change, the average surface temperature of the Earth will increase by Intergovernmental Panel on 
between 1.4°C and 5.8°C between 1990 and 2100. This compares with an increase Climate Change (2001) Climate 
of 0.6°C during the 20th century. Change 2001: Third Assessment 


Report. Cambridge University 
It is expected that this will be accompanied by changes in precipitation, retreat of | Press, Cambridge. 
glaciers and ice sheets in some parts of the world, and rising sea levels. 


Unsurprisingly, climate change on this scale has momentous implications for 
human health. Direct effects on health include, for example, changes in 
temperature-related mortality, and loss of life in storms and floods. Indirect 
effects are likely to be even more important, for example through reduced food 
and water availability and quality, and increased spread of some infectious 
diseases such as malaria. Lower-income populations, predominantly within 
tropical and subtropical countries, are likely to be the most affected. 


In this section, the direct impact of temperature on human mortality is explored. 
Thus, throughout the section, the outcome of interest is death, and the exposure 
is temperature. The relationship between temperature and mortality is 
complicated, as shown in Activity 1.1. 





Activity 1.1 The relationship between temperature and mortality 


The human body acts as an efficient heat regulator, keeping the internal 
temperature constant in different environments. As the ambient temperature 
reaches extremes of hot or cold, the body has to work harder to maintain the 
internal temperature, resulting in thermal stress. This stress in turn has adverse 
consequences for human health. 





A study was undertaken in Taiwan to estimate how the risk of death due to Pan, W.H., Li, L.A. and Tsai, 
various causes varies with ambient temperature. Mortality data were collected for M.J. (1995) Temperature 

the years 1981-1991 and related to ambient temperatures. For each temperature t ©Xtremes on mortality a 
(in °C), the risk of death was calculated for all days for which the ambient aL eee a te 
temperature was t. The minimum mortality temperature is the temperature at Chinese. Lancet, 345, 353-355. 
which the risk of death is lowest. Odds ratios were then calculated for each value 


of t, relative to the minimum mortality temperature. 


Figure 1.1 shows a line plot of the odds ratios by temperature, for death from 
cerebral infarction in persons aged 65 years or over. 


Odds ratio 
2 


8 12 16 20 24 28 a 
Temperature (°C) 


Figure 1.1 Odds ratio for death from cerebral infarction by temperature (°C) 


(a) Describe the dose-response relationship between ambient temperature and 
mortality in persons aged 65+ years due to cerebral infarction. 
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(b) Use Figure 1.1 to estimate the minimum mortality temperature in this 
population. 


(c) Briefly discuss the implications of Figure 1.1 for mortality if global 
temperatures increase. 


A dose-response curve that goes down then up as in Figure 1.1 is called U-shaped 
(even though in this case the U is rather asymmetric). U-shaped dose-response 
relationships between ambient temperature and mortality have been observed in 
different age groups, in different locations, and for different causes of death. 
However, the value of the minimum mortality temperature varies from place to 
place: it tends to be lower in colder climates than in warmer climates. For 
example, in the Netherlands it is 16.5°C, which is much lower than the value 29°C 
you found in Activity 1.1. The variation in minimum mortality temperatures is 
due to the fact that people adapt physiologically and culturally to their 
environment. 


Although the minimum mortality temperature varies from place to place, the 
strength of association between mortality and increases in temperature above the 
minimum mortality temperature can be compared for different locations. This is 
the topic of Activity 1.2. 


A comparative study was undertaken of the association between mortality and an 
increase of 1°C above the minimum mortality temperature in different locations. 
In one analysis, results obtained from five studies of mortality in persons aged 
65+ years from cardiovascular diseases (mainly heart attacks and strokes) were 
compared. Studies 1 and 2 used the odds ratio, while studies 3, 4 and 5 used the 
relative risk. In this activity, the term ‘effect estimate’ relates either to the 
relative risk or to the odds ratio, as appropriate. 





As part of the comparison, a meta-analysis was undertaken. The forest plot is 
shown in Figure 1.2. 








These data were adapted from 
Martens, W.J.M. (1998) 
Climate change, thermal stress 
and mortality changes. Social 
Science and Medicine, 46, 
331-344. 
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Figure 1.2 Forest plot for the association between mortality due to cardiovascular diseases and a 1°C rise above 


the minimum mortality temperature 
(a) In what circumstances is it appropriate to compare odds ratios and relative 
risks in this way? 


(b) Briefly compare the results from the five studies and the pooled values. 
Which study contributed most to the meta-analysis? 
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(c) A test for homogeneity similar to Tarone’s test (and with the same degrees of 
freedom) gives the value 20.27 for the test statistic. Obtain an approximate 
p value for the test, and interpret it. 





(d) Briefly discuss your results. Is it sensible to pool the effect estimates from 
these five studies? 


A U-shaped relationship with temperature has also been found for mortality due 
to respiratory diseases, and for total mortality. The estimated association between 
temperature and mortality is influenced to some degree by people taking action to 
reduce the impact of high or low temperatures — for example, by altering activity 
levels and clothing. The very young and the very old are less able to mitigate the 
effect of extreme temperatures in this way. In particular, there has been much 
interest in a possible association between sudden infant death syndrome (SIDS) 
and the thermal environment. This is the topic of Activity 1.3. 





A case-control study was undertaken in Tasmania of the association between 
sudden infant deaths and the thermal environment. Cases were infants who died 
at home between birth and one year of age and were later diagnosed to have died 
of SIDS. Controls were infants who survived to one year of age. The room 
temperature and the total thermal resistance of the baby’s clothing and bedding 
were measured — at the scene of death for the cases, and at the scene of last sleep 
for the controls. Thermal resistance is measured in togs: the higher the tog value, 
the greater is the thermal resistance. ‘These measurements were used to calculate 
the excess thermal insulation (in togs), over and above the recommended value for 
each room temperature. 


Cases and controls were classified as unexposed if their excess thermal insulation 
was within the range —4 to +4 togs, and exposed if their excess thermal 
insulation was greater than 4 togs. The data are in ‘Table 1.1. 


Table 1.1 Sudden infant deaths and excess thermal insulation 


Exposure Cases Controls 
Exposed: excess thermal insulation 9 4 
greater than or equal to 4 togs 

Unexposed: excess thermal 18 50 
insulation between —4 and 4 togs 

Total Zi 54 


(a) Estimate the odds ratio for the association between exposure and SIDS. 
(b) Obtain a 95% confidence interval for the odds ratio. 


(c) It is known that putting a baby down to sleep on its front increases the risk 
of SIDS. The authors comment that ‘sleeping position was found not to be a 
confounder of the relation between excess thermal insulation and SIDS’. 
Explain what this means, and how you might investigate the presence of 
confounding by sleeping position. 








(d) Summarize your results. 


The data from the case-control study described in Activity 1.3 can also be used to 
study the dose-response relationship between excess thermal insulation and SIDS. 





Ponsonby, A.L., Dwyer, T., 
Gibbons, L.E., Cochrane, J.A., 
Jones, M.E. and McCall, M.J. 
(1992) Thermal environment 
and sudden infant death 
syndrome: case-control study. 
British Medical Journal, 304, 
217-282. 


The tog value of a fabric is 
defined as 10 times the 
temperature difference in °C 
between its two faces when the 
heat flow is 1 watt /m?. 
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Activity 1.4 Thermal insulation and SIDS: dose-response analysis 


Excess thermal insulation can be quantified using tog values (see Activity 1.3). 
The value zero corresponds to thermal balance. Positive values correspond to high 
levels of thermal insulation, whereas negative values correspond to low levels of 
thermal insulation. Table 1.2 shows the numbers of cases and controls by dose, for 
excess tog values greater than or equal to —2. 


Table 1.2 Excess thermal insulation (togs) 
and SIDS, by dose 


Excess tog value Cases Controls 
> 6 5 1 
4 to 5.9 4 3 
2 to 3.9 5 {l 
0 to 1.9 í 13 
—2 to —0.1 4 21 


(a) Obtain the dose-specific odds ratios, relative to the —2 to —0.1 tog category, 
and comment briefly on your results. 





(b) The chi-squared test for no linear trend between level of excess thermal 
insulation and SIDS gives the value 10.026 for the test statistic. Obtain an 
approximate p value for the test, and interpret it. 


(c) There were three cases and five controls with excess tog values below —2. 
Obtain the dose-specific odds ratio for tog values below —2, relative to the 
—2 to —0.1 tog category. What does this suggest? 


(d) Briefly summarize your conclusions. 


To some extent, human physiology can adapt to gradual changes in temperature, 
as shown by the variation in the minimum mortality temperature between 
different locations on Earth. Furthermore, in temperate zones, the reduction in 
numbers of cold-related deaths resulting from a gradual increase in mean 
temperatures is likely to be greater than the increase in numbers of heat-related 
deaths. This suggests that the gradual increase in mean temperatures that is 
forecasted is unlikely, on its own, to lead to a substantial increase in mortality. 
However, it is also forecasted that there will be an increase in variability of 
temperatures, with more frequent heatwaves in particular. Human physiology is 
much less well adapted to coping with rapid changes in temperature than with 
gradual changes. 








Between 12 July 1995 and 16 July 1995, a heatwave occurred in Chicago, with 
temperatures reaching unprecedented highs; this was combined with very high 
relative humidity. A case-control study was undertaken to identify some of the 
risk factors associated with death. Cases were persons aged 24 years or more who 
died in Chicago between 14 July and 17 July, with a cause mentioned on the 
death certificate that was possibly heat-related. For each case, a matched control 
was selected, of the same age and living in the same neighbourhood. 





The risk factor of interest is participation in group activities involving social 
interactions. The data from this 1-1 matched case-control study are in Table 1.3. 


Semenza, J.C., Rubin, C.H., 
Falter, K.H. et al. (1996) 
Heat-related deaths during the 
July 1995 heat wave in Chicago. 
New England Journal of 
Medicine, 335, 84-90. 
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Table 1.3 Heat-related death and participation in group activities 





Controls 
Participated Did not participate 


Participated (1 63 


Cases _ 
Did not participate 90 14 





(a) Obtain the Mantel-Haenszel odds ratio for the association between dying of a 
heat-related disease and participation in group activities, and its 95% 
confidence interval. 


(b) Test the null hypothesis of no association between participating in group 
activities and dying of a heat-related disease, using McNemar’s test. 


(c) Summarize your findings. 


The study described in Activity 1.5 found that other factors associated with 
increased social contact, such as access to transport and even ownership of a pet, 
were also associated with a reduction in heat-related mortality. Owning or having 
access to a working air-conditioner was associated with a large reduction in 


heat-related deaths. However, use on a large scale of 
Se air-conditioning in the home is 
In concluding Section 1, it is important to emphasize that only the direct impact not sat hout adverse 


of increased temperatures on mortality has been considered. The impact of a environmental impact. 
gradual increase in mean temperature is likely to be relatively benign, at least in 

temperate zones. An increase in the frequency of heatwaves is likely to increase 

overall mortality. However, the indirect health impacts of global climate change 

— for example, concerning the availability of water and food — are likely to be 

much more serious, and are more difficult to quantify accurately. 














Summary of Section 1 


The dose-response relationship between ambient temperature and mortality is 
U-shaped. The minimum mortality temperature varies between regions of the 
world, owing to adaptation of human populations to their environment. A 
moderate increase in mean temperatures will result in a decrease in cold-related 
deaths and an increase in heat-related deaths, and a gradual shift in minimum 
mortality temperatures. However, human beings may not cope well with increased 
frequency of heatwaves. 








2 The changing temperature of the Earth 


The key observation behind current concerns about global climate change is that 
the average surface air temperature of the Earth is increasing, and increasing 
rapidly when compared to earlier periods. There are many ways in which this 
overall increase is manifested. These include local increases in average 
temperature (as for the central England temperatures series discussed in Book 2), 
and effects such as reductions in the extent of permafrost and sea ice in arctic 
regions, a rise in sea levels (due to both heat expansion and melting of sea ice), 
and changes in weather patterns. 











In this section, time series data on the overall change in the average surface air 
temperature of the Earth are considered. Calculating this average is not simple, 
as temperatures vary from place to place, and also vary with altitude, wind 
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conditions, type of land cover, and many other factors. For example, it is not even 
clear what is meant by ‘surface air temperature’ in a rain forest, since the 
temperatures at ground level and in the tree canopy will typically differ. There 
are further problems with how to measure surface air temperatures over the sea. 


However, while absolute temperatures are difficult to define and measure, 
temperature changes are much more robust to local characteristics. In climate 
science, the difference in temperature relative to an average value obtained over a 
defined baseline period is called an anomaly. Thus time series data on global 
temperature change are obtained by averaging data on temperature anomalies 
relative to the same baseline period from very many locations. In practice, various 
adjustments to the data are made prior to averaging, such as estimating and 
removing the ‘urban effect’ due to the growth of cities. (The growth of cities is 
estimated to have contributed no more than 0.05°C to the warming of the planet 
in the 20th century.) 





One of several available time series on global temperature anomalies is described 
in Activity 2.1. 


Activity 2.1 Global temperature anomalies, 1880-2006 


Monthly global temperature anomalies, relative to monthly averages over the 

30 years 1951-1980, were obtained from land-based meteorological stations for the 
period January 1880 to June 2006. Sea surface temperatures are not included in 
these data. 


(a) Figure 2.1 shows a seasonal plot of the monthly anomalies, for one year in 
each decade, namely 1880, 1890,..., 2000. 


Temperature anomaly (°C) 





8 
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec 
Month 


Figure 2.1 Seasonal plot of monthly global temperature anomalies 


Briefly describe any patterns in this plot. Is there any evidence of annual 
seasonality? Explain why or why not. 


10 


The term anomaly does not 
imply that the temperature 
difference to which it refers is 
anomalous. 


Hansen, J., Ruedy, R., 

Glascoe, J. and Sato, M. (1999) 
GISS analysis of surface 
temperature change. Journal of 
Geophysical Research, 104, 
30997-31022. 

These data were obtained from 
http: //data.giss.nasa.gov/gistemp 
in August 2006. 
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(b) A time plot of the data is shown in Figure 2.2. 


Anomaly (°C) 
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Figure 2.2 Time plot of global temperature anomalies 


Briefly describe the main features of this time plot. Is the time series 
stationary? 


(c) Figure 2.3 shows three moving averages of the data in Figure 2.2. 


Moving avcrage Moving avcrage Moving average 





0.5 
i a 
eee 
—(.5 J a 
Jan Jan Jan Jan Jan Jan Jan Jan Jan Jan Jan Jan Jan Jan Jan 
1880 1910 1940 1970 2000 1880 1910 1940 1970 2000 1880 1910 1940 1970 2000 
(a) Mouth (b) Mouth (c) Mouth 


Figure 2.3 ‘Three moving averages for the global temperature anomalies 


The orders of the moving averages are 61, 121 and 361 (months), 
corresponding to roughly 5, 10 and 30 years. 


Which moving average is which? In your view, which moving average gives 
the best estimate of the underlying trend? Explain your answer, and describe 
the trend. 





The time plots in Figures 2.2 and 2.3 suggest that the global average temperature 
increased by about 0.7°C between 1900 and 2000. An ARIMA model may be used 
to predict future increases. The selection of an ARIMA model for the global 
anomalies data is the subject of Activity 2.2. 


11 
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Activity 2.2 An ARIMA model for global temperature anomalies 


In part (b) of Activity 2.1, you saw that the variability in the data may be greater 
at the beginning of the monthly time series of global temperature anomalies than 
at the end, so the data for 1880 to 1899 will not be used. 


(a) In part (b) of Activity 2.1, you found that the time series of anomalies is not 
stationary. Time plots of the series of first and second differences of the 
monthly anomalies for the period January 1900 to June 2006 are shown in 
Figure 2.4. 


First differcnee Second difference 





Feb Feb Feb March March March 
1900 1950 2000 1900 1950 2000 
(a) Mouth (b) Mouth 


Figure 2.4 Global temperature anomalies: (a) first differences (b) second differences 


Identify the order of differencing required to obtain stationarity in mean. Is 
the differenced series stationary? Explain your answer. 





(b) Figure 2.5 shows the correlogram and partial correlogram for the first 
differences of the anomalies data. 


Autocorrclation Partial autocorrclation 
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(a) (b) 


Figure 2.5 First differences: (a) correlogram (b) partial correlogram 


(i) The length of the data series is n = 1277. Calculate the positions of the 
significance bounds in Figure 2.5. 


(ii) Suggest one or more possible ARIMA models for the anomalies data. In 
each case, explain your choice. 
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Activity 2.3 Predicting future temperature rises 


(a) An ARIMA(0,1,1) model was fitted to the anomalies data. Figure 2.6 shows 
the correlogram of the 1-step ahead forecast errors. 


Autocorrelation 
1 
0 — E E e a OEE ee 





12345 6 7 8 9 1011 1213 14 15 16 17 18 19 20 
Lag 


Figure 2.6 Correlogram of 1-step ahead forecast errors 


Interpret this correlogram. The p value of the Ljung—Box test of zero 
correlation at lags 1 to 20 was 0.001. Is the ARIMA(O, 1,1) model adequate? 
Explain your answer. 


(b) An ARIMA(1, 1,2) model was fitted to the data and was found to fit the 
data adequately. The parameter values for the fitted model were as follows: 


i ~ 0.001, 8, ~ 0.840, 
0, ~ 1.458, 05 ~—0.468, G~0.151. 
(i) Write down the model formula for the series of first differences. 


Figure 2.7 shows a multiple time plot of the observed anomalies to June 2006, 
the forecasts to December 2050 (obtained using the ARIMA(1, 1,2) model), 
and 95% prediction limits. 


Anomaly (°C) 





observed anomaly 
forecast 





— 95% prediction limits 








Mouth 


Figure 2.7 Observed anomalies, forecasts and 95% prediction limits 





(ii) Briefly describe the predictions, and the assumptions upon which they 
are based. 
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In Activity 2.3, an ARIMA model was used to extrapolate current trends into the 
future. ARIMA models do not incorporate any information about the mechanisms 
that give rise to variation in temperatures. In the case of global warming, much is 
known about the factors that influence the temperature of the Earth. These 
include solar activity, volcanic eruptions, increasing urbanization, and the 
emission of greenhouse gases. 





Greenhouse gases, which include carbon dioxide and methane, are gases that trap 
the heat of the Sun, like the glass panes of a greenhouse, thus resulting in an 
increase in air temperature. Changes in levels of greenhouse gases are believed to 
be the main cause of global warming. 


Activity 2.4 The growth in greenhouse gas emissions 


Carbon dioxide (CO2) is perhaps the single most important greenhouse gas. It is 

produced naturally, but also by the combustion of fossil fuels such as coal and oil. 

Emissions of CO have increased markedly since the beginning of the industrial 

era. Figure 2.8 shows the time plot of the logarithms of total annual emissions 

from fossil fuel burning, cement manufacture and gas flaring, in millions of tons of 

carbon. These data were obtained in 
August 2006 from the Carbon 
Dioxide Information Analysis 
Center website 

http: //cdiac.ornl.gov. 





log(CO:) 





5 
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Year 


Figure 2.8 Logarithms of total annual emissions of CO2, 1880-2003 
(a) Briefly describe the time plot in Figure 2.8. 


(b) Forecasts of total COə emissions are to be obtained using exponential 
smoothing. Explain why Holt’s method is appropriate, but simple 
exponential smoothing and the Holt—Winters method are not. 


(c) The first three values of the time series of logarithms are 5.46, 5.49 and 5.55. 
Obtain starting values for Holt’s method. 


(d) Holt’s method yields the optimal parameter values a = 1 and y = 0. 
Interpret these values. 


(e) Using Holt’s method, the forecasted value of the logarithm of the total 
emissions of CO» for 2050 is 10.3060. Calculate the forecasted total emissions 
of CO» for 2050, and briefly discuss the key assumption upon which this 
forecast is based. 


Climate scientists have a good understanding of how factors such as solar activity, 
volcanic eruptions, urbanization, and greenhouse gas emissions affect the 
mathematical equations they use to describe the atmosphere. From these 
equations, deductions can be made about how these factors influence the surface 
air temperature of the Earth. Thus scientists can use mathematical models to 
predict changes in the Earth’s temperature. Such predictions are the subject of 
Activity 2.5. 
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Activity 2.5 Natural and anthropogenic causes of global warming 


The factors affecting global warming may be categorized as natural and 
anthropogenic. Natural causes are those that would have occurred irrespective of 
human activity, such as changes in solar activity, and volcanic eruptions. 
Anthropogenic causes are those attributable to human activity, such as emissions 
of greenhouse gases from industrial and domestic sources. 


In order to evaluate the involvement of natural and anthropogenic factors in 
global warming, mathematical models are used to forecast global air 
temperatures, and the forecasts of the models are then compared with the 
observed data, in much the same way as for statistical models. 








Figure 2.9 shows time plots of observed temperature anomalies between 1860 and 
2000, relative to the first decades of the series, together with forecasts from three 
mathematical models. 
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Figure 2.9 Observed and forecasted temperature anomalies under contrasting models: (a) natural factors 
(b) anthropogenic factors (c) natural and anthropogenic factors. 


Figure 2.9(a) shows forecasts based on natural factors only, namely variation in Figure 2.9 was adapted from 

solar activity, and volcanic eruptions. Figure 2.9(b) shows forecasts based on Figure SPM-2 of the following 

anthropogenic factors only, primarily emissions of greenhouse gases. Figure 2.9(c)  TePort: Climate Change 2001: 
Summary for Policymakers, 

shows forecasts based on both natural and anthropogenic factors. Copyright © 2001, 

(a) Use Figure 2.9 to comment on the accuracy of the forecasts obtained using Intergovernmental Panel on 


Climate Change, Cambridge 


the three models, and hence identify the best model. a 
University Press. 


(b) What does the model that you identified in part (a) suggest about the likely 
role of anthropogenic factors in global warming? 
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Summary of Section 2 


Changes in the mean surface air temperature of the Earth may be measured by 
averaging temperature anomalies from many locations. One data set on global 
temperature anomalies suggests that the mean surface air temperature has 
increased by about 0.7°C between 1900 and 2000. Emissions of greenhouse gases 
have grown steeply over the same period. Comparisons of observations and 
forecasts from mathematical models suggest that both natural and anthropogenic 
factors affect global temperatures, but that anthropogenic factors — and, 
specifically, the increase in greenhouse gas emissions — have contributed strongly 
to the increase in mean temperatures since the 1970s. 


3 Predicting weather patterns 





It is predicted that, as well as the increase in average temperatures, weather 
patterns on the Earth will become more variable in future. There is already 
evidence that such variability is increasing, as shown by more frequent floods, 
droughts, heatwaves and cold spells in many parts of the world. In addition, the 
impact of climate change is likely to vary between regions. 





To understand the impact of global warming better, it is necessary to understand 
the mechanisms that drive the weather in different parts of the world. Weather 
data typically comprise measurements on several variables — for example, 
temperature, atmospheric pressure, humidity, precipitation, wind strength and 
wind direction — taken at many locations, at many time points and at several 
levels of the atmosphere. The techniques of multivariate analysis may be used to 
identify the dominant patterns in such data. 





Oceans are particularly important in shaping weather patterns. In this section, 
weather patterns related to sea-level atmospheric pressure are considered. 


Activity 3.1 Sea-level atmospheric pressure over the North Atlantic 


The weather in western Europe, and in the UK in particular, is greatly influenced 
by proximity to the Atlantic Ocean. Data were obtained for the years 1865 to 
2001 on atmospheric pressure at sea level at three locations — Reykjavik in 
Iceland, Ponta Delgada in the Azores, and Gibraltar (close to the southern tip of 
Spain). These locations are shown in Figure 3.1. 


The mean atmospheric pressures in the winter months (January to March) and in 
the summer months (June to August) were obtained for each year at each of the 
three locations. Hence there are six variables (two measurements at each of three 
locations). The observations are the values obtained in different years. Figure 3.2 
shows a matrix scatterplot of the data. 





(a) Briefly describe the relationships, if any, between the six variables. Are the 
relationships between measurements taken at different locations stronger in 
winter or in summer? Explain your answer. 


(b) The diagonal and lower triangle of the correlation matrix for the data are 
given below. 


Iceland, winter 1 

Azores, winter —0.712 1 

Gibraltar, winter —0.646 0.642 1 

Iceland, summer 0.152 —0.106 —0.026 1 

Azores, summer 0.057 0.023 —0.010 —0.357 1 


Gibraltar, summer —0.098 0.084 0.081 0.170 0.192 1 
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The data and information in 
this section were obtained in 
August 2006 from the website of 
the Climatic Research Unit 
http://www.cru.uea.ac.uk. 


Section 3 Predicting weather patterns 


Comment on the relative strength of the relationships between measurements 


made at the same location at different times of year, and at different 
locations at the same time of year. 
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Figure 3.1 The North Atlantic 
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Figure 3.2 Matrix scatterplot of sea-level pressure data 
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Activity 3.2 Pressure patterns over the North Atlantic 


In Activity 3.1, you found that atmospheric pressures in Iceland are negatively 
related to those in Gibraltar and the Azores. Thus when the pressure is high in 
Iceland, it tends to be low in Gibraltar and the Azores, particularly in the winter 
months. 


Principal component analysis (PCA) is commonly used to describe the patterns in 
sea-level atmospheric pressure data obtained from several locations. In this 
activity, PCA is applied to the data on pressures. 





(a) The atmospheric pressures were measured in millibars. The mean and 
variance of the six pressure variables are given in Table 3.1. 


Table 3.1 Mean and variance of sea-level pressures 


Location and season Mean Variance 
Iceland, winter 1001.71 33.08 
Azores, winter 1020.08 14.72 
Gibraltar, winter 1020.63 5.55 
Iceland, summer 1009.79 4.88 
Azores, summer 1024.14 1.75 


Gibraltar, summer 1016.36 0.37 


(i) Briefly discuss these results. 


(ii) Why might it be advisable to base a principal component analysis on 
standardized data? 


(iii) Explain why at most six principal components can be obtained from 
these data. 


(b) Figure 3.3 shows the scree plot obtained in the principal component analysis. 


Variance 





Il 2 3 4 J 6 
Component number 


Figure 3.3 Scree plot for the sea-level pressure data 





(i) Based on the scree plot, how many principal components should be 
retained? 


(ii) How many principal components should be retained according to 
Kaiser’s criterion? 





(c) It is decided to retain three principal components. Their variances are given Table 3.2 Variances of 
in Table 3.2. principal components 
Calculate the percentage variance explained by each of these components, oe l 

; , : component Variance 
and the cumulative percentage variance explained by all three components. — neo 

(d) Briefly summarize the results of the analysis so far. 2 1.351 

3 1.131 
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In Activity 3.2, you saw that the atmospheric pressure data can be approximated 
by a data set of lower dimension. The first two principal components represent an 
important feature which affects the weather in western Europe. This feature is 
called the North Atlantic Oscillation. 


Activity 3.3 The North Atlantic Oscillation 


In Activity 3.2, you found that the six-dimensional data set on winter and 
summer pressures can be adequately approximated by a three-dimensional data 
set based on the first three principal components of the standardized data. The 
loadings of these principal components are given in Table 3.3. 


Table 3.3 Principal component loadings 


Principal component 


Location and season I 2 3 

Iceland, winter —0.580 0.049 0.063 
Azores, winter 0.577 —0.025 —0.027 
Gibraltar, winter 0.553 —0.101 —0.002 
Iceland, summer —0.114 —0.687 0.361 
Azores, summer 0.025 0.718 0.319 


Gibraltar, summer 0.099 0.017 0.874 
Component variance 2.368 1:391 1131 


(a) Interpret the three principal components. 


(b) Figure 3.4 shows a scatterplot of the first two principal components. 


Second principal component 





First principal component 


Figure 3.4 Scatterplot of the first two principal components 


(i) Explain why you would not expect the scatterplot to show any evidence 
of a linear relationship between the two components. 


(ii) Briefly explain what features distinguish the points corresponding to 
1902 and 1990. 


(c) Figure 3.5 (overleaf) shows a scatterplot of the mean January to March 
temperature in central England (in °C) against the value of the first principal 
component, for the years 1865 to 2001. 
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First principal component 


Figure 3.5 Scatterplot of mean January to March central England temperature 
against first principal component 


(i) Comment on the relationship between the mean January to March 
temperature in central England and the value of the first principal 
component. To which year does the point in the top right-hand corner 
of the scatterplot correspond? 


(ii) The correlation between the first principal component and the mean 
January to March temperature is 0.723. Interpret this correlation. 


In Activity 3.3, you saw that the first principal component represents the 
difference between the (standardized) mean sea-level pressures measured in 
winter, in the Azores and Gibraltar on the one hand, and Iceland on the other. 
The second principal component is also a south minus north difference between 
mean pressures, but in summer. You also found a strong positive association 
between the first principal component and the mean winter temperature in central 
England. 











Climate scientists use such pressure differences to monitor a key feature of weather 
patterns known as the North Atlantic Oscillation (NAO). This determines the 
frequency and path of storms over the North Atlantic, and hence the weather 
(including temperature and rainfall) in different parts of western Europe. For 
example, in high NAO winters (corresponding to large south minus north pressure 
differences over the North Atlantic, and hence to large positive values of the first 
principal component), storms bring warm wet air over north-western Europe, 
resulting in abnormally warm and wet winters. Conversely, during low NAO 
winters (corresponding to low south minus north pressure differences, and hence 
to low values of the first principal component), the storms move further south and 
the weather in north-western Europe tends to be cold and dry. 








Climate scientists monitor the NAO using a simple index based on the 
standardized sea-level pressure difference in winter months between Gibraltar and 
Iceland, or between the Azores and Iceland. How effective this simple index is at 
capturing temperature variations in central England is the subject of Activity 3.4. 





Activity 3.4 The NAO index 


Mean January to March temperatures in central England for the years 1865-2001 
were grouped into three categories based on the 0.33-quantile and the 
0.67-quantile: cold, average and warm. A discriminant analysis was undertaken 
using the group-standardized pressure measurements in the winter months 
January to March for Iceland (Z1), the Azores (Za) and Gibraltar (Za). 


(a) How many useful discriminant functions can be obtained? Explain your 
answer. 
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More elaborate indices, derived 
using principal components, are 
also commonly used to monitor 


the NAO. 


Section 3 Predicting weather patterns 


(b) The first discriminant function is 
D = 0.676Zg — 0.7722; — 0.346 Z4. 


Interpret this discriminant function. 





(c) The NAO index, represented by the standardized Gibraltar pressure minus 
the standardized Iceland pressure, has the characteristics shown in Table 3.4. 


Table 3.4 Group means and variances 
for the NAO index 


Group n Mean Variance 
Cold 46 —1.6083 1.9031 
Average 46 0.0446 1.4736 


Warm 45 1.5984 1.4033 


Obtain the between-groups variance and the within-groups variance, and 
hence calculate the separation achieved by the NAO index. 


(d) The separation achieved by the first discriminant function D is 1.170, and the 
total separation achievable is 1.171. Calculate the percentage separation 
achieved by D and by the NAO index. How effective is the NAO index in 
separating the groups, compared to D? 





The NAO index provides a useful indicator of major weather patterns in western 
Europe. A similar feature in the southern hemisphere is the subject of 
Activity 3.5. 


Activity 3.5 El Niño events and the Southern Oscillation 


Fishermen off the coast of Peru have long noticed that in some years the waters of 
the Pacific Ocean are uncommonly warm, particularly over the Christmas period. 
They call this phenomenon ‘El Nino’, which in Spanish means ‘the boy’ — a 
reference to the newborn Christ. Unusually high ocean temperature in the East 
Pacific is just one of many attributes of what have since been called El Nino 
events, which are also associated with abnormal fluctuations in ocean currents and 
atmospheric conditions. El Nino events have consequences for weather patterns in 
areas beyond the Pacific Ocean, including India, Africa and North America. 


The occurrence of El Nino events is closely associated with variations in pressure 
in different parts of the Pacific Ocean. This variation is known as the Southern 
Oscillation. A map of the equatorial Pacific Ocean is shown in Figure 3.6. 
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Figure 3.6 The equatorial Pacific Ocean 


The first discriminant function 
is given in part (b). 


The opposite phenomenon is 
called ‘La Nina’ — the girl. 
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Figure 3.7(a) shows a scatterplot of annual mean sea-level pressures (in millibars) 
in two locations, Tahiti (an island in the centre of the Pacific) and Darwin (a 
coastal city in northern Australia), for the years 1866-2005. 


Figure 3.7(b) shows a scatterplot of the monthly mean pressures for the same 
period. 
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Figure 3.7 Scatterplots of atmospheric pressures (in millibars) in Darwin and Tahiti: (a) annual means 


(b) monthly means 


(a) Describe the relationship between the annual mean pressures in Darwin and 
in Tahiti, and the relationship between the monthly mean pressures in 
Darwin and in Tahiti. What do you notice? 


(b) Suggest an explanation for your finding in part (a). You may find it useful to 
think of the data in Figure 3.7(b) as comprising twelve groups, group 1 
containing the January points for the 140 years 1866-2005, group 2 
containing the February points, and so on. Then consider the relationships 
within-month and between-month (that is, within and between the monthly 


groups). 


(c) The between-month covariance matrix B and the within-month covariance 
matrix W for monthly pressures at Darwin (X,) and Tahiti (X2) are 


5.80 3.00 LIQ =0}32 
oe o 4 - E a) i 
Do these matrices support the explanation that you suggested in part (b)? 
Explain your reasoning. 





(d) The annual Southern Oscillation index is defined as the difference between 
the annual mean sea-level pressures in Tahiti and in Darwin, suitably 
standardized. Explain in general terms why this is a sensible one-dimensional 
approximation to the data in Figure 3.7(a). 


El Niño events typically occur once every four to seven years, and last about one 
or two years. There is evidence that the frequency and duration of El Niño events 
increased in the latter part of the 20th century, compared to earlier times. One 
way to study the impact of global warming on weather patterns is to monitor 
changes over time in indices such as the Southern Oscillation index and the North 
Atlantic Oscillation index. For example, there is some evidence that the value of 
the NAO index increased in the period 1960-1990. In order to better understand 
the implications of global warming for weather patterns across the world, it is 
important to gain a better understanding of how it affects variables such as the 
North Atlantic Oscillation index and the Southern Oscillation index. 
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The explanation is quite subtle. 
You may wish to look at the 
solution before trying part (c). 


Section 4 Changing sea levels 


Summary of Section 3 


The effect of global warming on weather patterns in different parts of the world is 
likely to be very variable, and is hard to predict. To gain some insight into the 
likely effects of global warming, it is necessary to understand what factors affect 
weather patterns. Important features such as the North Atlantic Oscillation, 
which affects the weather in western Europe, and the Southern Oscillation, which 
is related to the occurrence of El Nino events, are closely associated with sea-level 
pressures at relevant locations. These pressures may be analysed using 
multivariate statistical methods to identify simple indices, such as the North 
Atlantic Oscillation index or the Southern Oscillation index, on which predictions 
may be based. 


4 Changing sea levels 


One of the most striking consequences of global warming is the rise in mean sea The data and information used 
levels. This has potentially serious consequences for many coastal or island in this section were obtained 
nations, and potentially devastating implications for low-lying regions, such as the '"0™ the following report. 


Intergovernmental Panel on 


atolls of Micronesia and the Marshall Islands (see Figure 3.6). Climate Change 001) Climate 


Change 2001: Third Assessment 
Report. Cambridge University 
Press, Cambridge. 


However, evaluating the impact of global warming on sea levels is particularly 
difficult, because sea levels also vary for reasons unconnected with climate change. 
Thus expert knowledge is required to assess the contributions of the different 
factors to the observed changes, and a Bayesian approach is required to combine 
data and expert opinions. 





In this section, changes in sea levels during the course of the 20th century are 

investigated. The factors that contribute to such changes may be classified as 

natural or anthropogenic (that is, related to global warming resulting from human These terms were introduced in 
activity). Natural factors include, for example, changes due to very long-term Activity 2.5. 

movements of the continental masses following the end of the last ice age, and 

deposition of sediments on the ocean floors. Anthropogenic factors, resulting from 

the impact of human activity on global temperatures, include thermal expansion, 

and the melting of permafrost, ice sheets and glaciers. 


Activity 4.1 Summarizing expert knowledge 


Scientists use a range of techniques to estimate the contribution of different 
factors to changes in sea levels. Some of these factors contribute a rise in sea 
levels, others contribute a reduction. ‘These estimates, and the uncertainties 
underpinning them, may be used to obtain prior distributions that summarize the 
current knowledge of experts. Two of these priors are shown in Figure 4.1 
(overleaf). 
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Figure 4.1 Priors for mean annual rise in sea levels due to natural factors and to 
anthropogenic factors 





The prior densities represent the beliefs of experts about the mean annual rise in 
sea levels (in mm/yr) during the period 1910-1990. Negative values indicate a 
decrease in mean sea levels. 


(a) Contrast the degree of uncertainty represented by these two priors. Which is 
the stronger prior? Explain your answer. 








(b) Briefly describe the experts’ beliefs about the contributions of natural and 
anthropogenic factors to changing sea levels between 1910 and 1990. 


(c) Experts’ beliefs about the mean annual rise in sea levels from the combined 
effects of natural and anthropogenic factors may be summarized using the 
following values (in mm/yr): 


mode = 0.70, 0.25-quantile = 0.41, 0.75-quantile = 0.99. 


(i) | Give two reasons why a normal prior N (a,b) may be appropriate for 
representing these beliefs. 


(ii) Use the assessed values of the mode and quartiles to specify a normal 
prior representing experts’ beliefs about the mean annual rise in sea 
levels due to natural and anthropogenic factors. 


The experts’ beliefs summarize experts’ best estimates of the contribution of 
various factors to rises in sea levels. These estimates are based on mathematical 
models that describe how glaciers and ice caps melt, how continental masses 
move, how the depth of permafrost in Arctic regions varies, how the volume of the 
oceans expands as the Earth surface warms up, and so on. However, the combined 
effect of all these factors may be measured directly, using observations on sea 
levels collected during the years 1910 to 1990. 


Activity 4.2 Combining sources of information 


Accurate measurements of mean changes in sea level are difficult to obtain. 
Measurements are usually made close to coastlines — for example, with tidal 
gauges. Such measurements are often affected by coastal features. Very few 
measurements of long-term changes in sea level from mid-ocean locations are 
available. Furthermore, measurements made from coastal locations measure only 
the change in sea level relative to the land, and must therefore be adjusted for 
vertical land movement. 








Thus direct measurements of changes in sea level may be biased, and are very 
variable. Nevertheless, they complement usefully the indirect calculations made 
by experts. In this activity, you will combine experts’ beliefs about the mean rise 
in sea levels with measurement data to obtain the posterior distribution for the 
mean annual rise in sea levels between 1910 and 1990. 
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(a) Measurements of the mean annual change in sea level may be regarded as 
observations from a normal distribution N (u, o°), with a? known. 


(i) Write down the posterior distribution for u, when the prior for py is the 
conjugate prior N (a,b). 


(ii) Show that the posterior mean and variance can be written in terms of 
a, b, and o?/n, where T is the sample mean of the measurements and 
n is the sample size, as follows: 


2 2 ha 
posterior mean = (o /n) x a+ OF 


(jJn) +b ” 

b 

posterior variance = CLA 

(o? /n) +b 

(b) The measurements yield the following summary statistics: 
z œ~ 1.50, — ~ 0.065. 

n 
2 


Assuming that g? is equal to s*, obtain the mean and variance of the 
posterior for u, using the prior that you obtained in Activity 4.1. 


The posterior distribution for the mean annual rise in sea levels summarizes all 
the available information from experts and data about the mean annual rise in sea 
levels during the period 1910-1990. In Activity 4.3, you are asked to consider the 
relationship between the posterior, the prior and the data. 


Activity 4.3 Summarizing the mean annual rise in sea levels 


Figure 4.2 shows the prior, the likelihood, and the posterior for u, the mean 
annual rise in sea levels during the period 1910—1990. 


--- likelihood 
—— posterior \ 
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Figure 4.2 Prior, likelihood and posterior for mean annual sea-level rise 


(a) Comment briefly on the location of the posterior in relation to the prior and 
the likelihood. 


(b) Explain how the measurement data have altered the prior uncertainty 
about u. 


(c) Comment on the posterior probability that there was a mean annual drop in 
sea levels during the period 1910-1990. 


(d) A 95% HPD credible interval for u is (0.86, 1.72). 
(i) Interpret this interval. 


(ii) Would the 95% equal-tailed credible interval for u be different? Explain 
your answer. 
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In Activities 4.2 and 4.3, a Bayesian analysis was undertaken to estimate the mean 
annual rise in sea levels during the period 1910-1990. This rise comprises two 
components: a rise due to natural factors, and a rise due to anthropogenic factors. 


In order to assess the contribution of global warming to changing sea levels, the 
magnitude of the rise due to anthropogenic factors is of greater interest. 
Unfortunately, measurements of changes in sea levels cannot differentiate between 
natural and anthropogenic causes. 








To see this, let z denote the sample mean of the measured mean annual rise in sea 
levels. This may be regarded as an observation from a normal distribution with 
mean u = H, + u, and variance o7/n, where u, is the mean annual rise due to 
natural factors, and u, is the mean annual rise due to anthropogenic factors. Thus 


— g? 
KON (ma Hio E) 


The parameter of interest is u,, which represents the mean annual rise due to 
anthropogenic factors. This cannot be estimated from measurements alone: 
inferences can be made only about u = u, + u,- In a classical (that is, 
non-Bayesian) framework, nothing more can be done. However, in a Bayesian 
framework, experts’ opinions can be used to evaluate the contribution of 
anthropogenic factors, and the data can be used to update this prior information 
about u„ and u,. A Bayesian model for u„ and u, is the subject of Activity 4.4. 





The following WinBUGS code specifies a model for the sample mean annual rise 
in sea levels, in terms of the parameters pu, and p,- 


model 
{ 
x ~ dnorm(mu, prec) 
prec <- 1/0.065 
mu <- mun + mua 
mun ~ dnorm(-0.075,precn) 
mua ~ dnorm(0.775,preca) 
precn <- 1/0.162 
preca <= 1/0.023 
} 


In this model, the variable (or node) x represents the sample mean of the 
measurements on the annual rise in sea levels. As described in Activity 4.2(b), 
there is a single observation on the sample mean, namely 1.50mm/yr. The 
variance of the sample mean is assumed known and equal to s?/n = 0.065. 


(a) The WinBUGS code involves seven nodes. Classify each node as either 
stochastic, deterministic or constant. 


(b) The nodes mun and mua represent the quantities u„ and u,, respectively. 
Write down the distributions for these parameters specified in the model 
definition. 





(c) Hence describe the specification of the model for the sample mean annual rise 
in sea levels. 
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These distributions were plotted 
in Figure 4.1. 


Section 4 Changing sea levels 


In the model defined in Activity 4.4, the data (namely, the sample mean annual 
rise in sea levels) provide information on the parameters p,, and u, by restricting 
the value of their sum. However, the prior information still plays the key role in 
determining the values of these parameters. MCMC was used to sample values 
from the posterior distribution for this model. The results of the simulation are 
the subject of Activity 4.5. 


Activity 4.5 Estimating the rise in sea levels attributable to 
anthropogenic factors 


WinBUGS was used to sample values from the model defined in Activity 4.4. 
Three chains of values were obtained, each of length 10000, using MCMC. ‘Trace 
plots of the values of u„ and u, for the three chains, superimposed on the same 
diagram, are shown in Figure 4.3. 





(a) Ha I 2500 5000 T500 10000 
Iteration 





(b) H, I 2500 5000 T500 10000 
Iteration 


Figure 4.3 Trace plots of simulated values: (a) u, (b) Hn 


(a) Explain whether or not a burn-in is required. 


(b) The prior mean and standard deviation for u„ and u,, and estimates of the 
posterior mean and standard deviation, are given in Table 4.1, together with 
the Monte Carlo standard error. 


Table 4.1 Prior and posterior 


Prior Posterior 
Prior standard Posterior standard 
Parameter mean deviation mean deviation MC error 
ia —0.075 0.4025 0.444 0.2385 0.0012 
[La 0.775 0.1517 0.847 0.1440 0.0001 


(i) The posterior mean and standard deviation are based on the values in 
all three chains, a total of 30000 simulated values. Use the 5% rule of 
thumb to decide whether sufficient values were simulated. 


(ii) Contrast the prior and posterior means and standard deviations, and 


hence comment on how the means and standard deviations of u„ and u, 


have been updated. For which parameter has the mean changed more? 
Suggest an explanation. 
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(c) The kernel density estimate of the posterior p.d.f. for u, is shown in 
Figure 4.4. 
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Figure 4.4 Estimated posterior density for pu, 


(i) Briefly describe the posterior distribution for p.. 


(ii) The estimated 0.025-quantile and 0.975-quantile of the posterior for u, 
are 0.563 and 1.128, respectively. Obtain an estimate of the 95% 
credible interval for u,. Is this an HPD interval or an equal-tailed 
interval? 





(iii) Summarize the likely contribution of anthropogenic factors to the mean 
annual rise in sea levels between 1910 and 1990. 


A Bayesian approach makes it possible to use both data and expert opinion on 
the rise in sea levels due to anthropogenic factors. The conclusion from this 
analysis is that global warming related to anthropogenic factors has resulted in a 
mean rise of about 0.85 mm per year over the period 1910—1990. 


Summary of Section 4 





There is evidence that sea levels rose in the 20th century. However, the 
contribution to this rise from anthropogenic factors related to global warming is 
more difficult to determine. Experts have developed models that predict the 
contribution of natural and anthropogenic factors to rising sea levels. The 
predictions from these models may be used as priors in a Bayesian analysis; the 
priors can be combined with the likelihood obtained from measurements of 
changing sea levels. The posterior distribution from this analysis suggests that 
anthropogenic factors contributed substantially to the mean annual rise in sea 
levels between 1910 and 1990. 
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5 Where next? 


In Sections 1 to 4, some of the methods discussed in Books 1 to 4 have been used 
to explore some statistical issues relating to global climate change. Of course, 
there is much more to climate science than could be covered here, and much more 
to modern statistics than could be covered in M249. In this section, some pointers 
are given to further areas of practical modern statistics. 





Medical statistics 


A central theme in medical statistics is the disentangling of the contributions of 
different variables, some of them confounders, in order to throw light on 
important associations, and hence, hopefully, on underlying causes of disease. The 
key method presented in M249 is stratification using the Mantel-Haenszel 
approach. However, there is a more general and more flexible (though perhaps 
less transparent) method called logistic regression. This is a particular instance of 
a general modelling approach known as generalized linear modelling, which itself 
can be further extended in various ways, for example to generalized mixed models. 








Another important area of medical statistics, which has not been covered in 
M249, is survival analysis. This is widely used in cohort studies of long duration, 
and when individuals are followed up for different lengths of time. In survival 
analysis, the emphasis is on a quantity known as the hazard function, which 
describes the rate at which events occur, rather than the probability of an event, 
over the duration of the study. In survival analysis, this rate is allowed to vary 
over time. Survival analysis can be used to analyse data on any type of life event 
including onset, relapse and progression of disease, as well as death. 


Time series 


The time series models that you have met in M249 provide a comprehensive 
introduction to time series modelling in what is known as the time domain, that 
is, to explain changes over time. Also important, but not covered in M249, is 
modelling in the frequency domain. Here the emphasis is on finding cyclical 
patterns in time series, and estimating the periods of such cycles, using a set of 
techniques known as spectral analysis. This is particularly useful in subject areas 
such as meteorology for understanding underlying regularities in long data series, 
beyond the well-known annual seasonal cycle. 


In the additive models discussed in M249, the variance is assumed to be constant. 
Many models have been proposed that allow for randomly changing variance, 
which in this context is called volatility. These include ARCH models 
(autoregressive conditional heteroscedasticity) and GARCH models (generalized 
ARCH), which are widely used in financial applications. 


Multivariate analysis 


The methods of multivariate analysis presented in M249 are all of a descriptive 
nature, in that they assume no underlying model for the data. Multivariate 
models have not been discussed in M249. The technique of factor analysis, for 
example, uses models based on the multivariate normal distribution. Factor 
analysis is superficially similar to principal component analysis, in that it seeks a 
reduction in dimension of the data. However, it assumes an underlying model, 
which PCA does not. Factor analysis is much used in areas such as psychology 
and marketing. 


A further area of multivariate analysis not covered in M249 is cluster analysis. In 
cluster analysis, the aim is to identify the groups, unlike in discrimination, in 
which the groups are known in advance. This is done using distance measures 
similar to the separation in discriminant analysis. Cluster analysis techniques are 
much used in taxonomy, in molecular and evolutionary biology and in ecology. 


Logistic regression and 


generalized linear models are 


described in M346 Linear 
Statistical Modelling. 


29 


Statistical perspectives on climate change 


Bayesian Statistics 


There has been a huge upsurge in the use of Bayesian methods since the early 
1990s, owing to the availability of fast computing techniques based on Markov 
chain Monte Carlo (MCMC). 


In M249, it has been possible to give only a flavour of the applications of modern 
Bayesian methods. In fact, Bayesian models applied using MCMC have opened 
up new fields of application for statistics — in particular, applications to problems 
requiring large numbers of parameters, problems involving missing data, and 
complicated applications in which the number of parameters in the model is not 
fixed, known as trans-dimensional problems. 


However, there is more to Bayesian analysis than complicated models and 
computer code. The theory of Bayesian inference has only been touched upon in 
M249. There is, for example, a Bayesian approach to testing hypotheses based on 
quantities known as Bayes’ factors, and a Bayesian theory of estimation based on 
quantities called loss functions. A further important topic of Bayesian statistics is 
decision analysis. In decision analysis, Bayesian methods are used to determine 
the best course of action by optimizing a quantity known as a utility function, 
which represents costs and benefits. 
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Markov chains, and many other 
stochastic models, are described 
in M343 Applications of 
Probability. 


Solutions to Activities 


Solution 1.1 


The odds ratio is defined in Subsection 1.2 of Book 1. 
Dose-response relationships are discussed in 
Subsection 8.2. 


(a) The odds ratio declines as the ambient 
temperature rises from 9°C to 29°C, then increases 
steeply as the temperature rises from 29°C to 32°C. 
Therefore the dose-response relationship is not linear. 





(b) The minimum mortality temperature is 29°C. 


(c) As average global temperatures increase, the 
number of deaths resulting from cold stress (that is, 
deaths resulting from temperatures lower than the 
minimum mortality temperature) will decline, and the 
number of deaths from heat stress (that is, deaths 
resulting from temperatures higher than the minimum 
mortality temperature) will increase. 


Solution 1.2 


The relationship between odds ratios and relative risks 
is discussed in Subsection 1.3 of Book 1. The 
interpretation of forest plots is described in 

Subsection 12.2. 


(a) For an uncommon disease, the relative risk and 
the odds ratio (for the same exposure) are 
approximately equal. Thus it is valid to compare odds 
ratios and relative risks obtained in different studies 
provided that deaths from cardiovascular diseases are 
uncommon in the populations studied. (However, 
cardiovascular diseases may not be uncommon in 
people over 65 years old.) 





(b) There is much variability between studies: the 
effect estimates were much higher for studies 2 and 5 
than for studies 1, 3 and 4. The confidence interval for 
the pooled estimate does not overlap with the 
confidence interval for study 2. Study 3 contributed 
most to the meta-analysis, as this study has the 
largest square. 


(c) Since there are five studies, the distribution of the 
test statistic under the null hypothesis of homogeneity 
is y7(4). From Table 3 of the Handbook, the 
0.999-quantile of y?(4) is 18.47. Since the observed 
value of the test statistic is greater than this, the 

p value is less than 0.001. Thus there is strong 
evidence of heterogeneity — that is, the strength of 
association between mortality and an increase of 1°C 
above the minimum mortality temperature varies 
between locations. 











(d) There is substantial heterogeneity between 
studies. For this reason, it is not sensible to pool 
results from the studies. 


Solution 1.3 


Confidence intervals for the odds ratio are discussed in 
Subsection 2.3 of Book 1. The analysis of case-control 
studies is described in Section 3. Confounding is 
defined in Subsection 6.4, and evidence for 
confounding is discussed in Subsection 7.1. 


(a) The estimated odds ratio is obtained using 
Formula (2.5) of Book 1: 








(b) For a 95% confidence interval, z = 1.96. The 
estimated standard error is given by Formula (2.7) of 
Book 1: 


Hence the 95% confidence limits are 
OR =ORx exp(—1.96 x a) 
~ 6.25 x exp(—1.96 x 0.6608) 
S (lL. 
OR” =0Rx exp(1.96 x a) 
~ 6.25 x exp(1.96 x 0.6608) 
D282, 


Thus a 95% confidence interval for the odds ratio is 
(defi, 22.82). 


(c) The association between SIDS and excess thermal 
insulation could be confounded by sleeping position if 
sleeping position were associated with excess thermal 
insulation (it is already known to be associated with 
SIDS). This could be investigated by undertaking an 
analysis stratified by sleeping position, and comparing 
the Mantel—Haenszel odds ratio with the crude odds 
ratio. If the two are similar then confounding is not a 
problem. 





(d) There is a strong positive association between 
SIDS and excess thermal insulation, odds ratio 6.25, 
95% confidence interval (1.71, 22.82). The association 
is not confounded by sleeping position. 
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Solution 1.4 

Dose-response analysis is discussed in Subsection 8.2 
of Book 1. 

(a) The dose-specific odds ratios are shown in 
Table $.1. 


Table S.1 ‘Thermal insulation and SIDS: 
dose-specific odds ratios 


Excess tog value Cases Controls Odds ratio 


> 6 5 1 26.25 
4 to 5.9 4 3 7.00 
2 to 3.9 5 11 2.39 
0 to 1.9 7 13 2.83 

~2 to —0.1 4 21 1.00 


The odds ratio tends to increase with excess tog value, 
indicating an increasing dose-response relationship 
between excess thermal insulation and SIDS. 


(b) The null distribution of the test statistic is y7(1). 
From Table 3 of the Handbook, the 0.995-quantile of 
x7(1) is 7.88 and the 0.999-quantile is 10.83, so 

0.001 < p < 0.005. Thus there is strong evidence of a 
linear dose-response relationship between the level of 
excess thermal insulation and SIDS, for excess thermal 
insulation in the range considered. 





(c) The dose-specific odds ratio is 

3x 21 

4x5 
This suggests that the risk of SIDS increases when the 
thermal insulation is more than 2 togs below the level 
at which thermal balance is achieved. 





= 3.15. 


(d) For excess thermal insulation in the range —2 to 
6+ togs, there is a positive dose-response relationship 
between excess thermal insulation and SIDS (test for 
no linear trend, p < 0.005). Insufficient thermal 
insulation (excess tog value less than —2) may also be 
positively associated with SIDS. This suggests that 
the dose-response relationship over the entire range 
may be U-shaped, with minimum mortality in the 
range —2 to 0 togs of excess insulation. 





Solution 1.5 


The analysis of 1-1 matched case-control studies is 
described in Subsection 7.2 of Book 1. 


(a) The Mantel-Haenszel odds ratio in a 1-1 matched 
case-control study is obtained using Formula (7.3) of 
Book 1: 

f 63 


OR =— = — =0.70. 
g 90 
For a 95% confidence interval, z = 1.96. The estimated 
standard error is calculated using Formula (7.4): 


C= a Pade 2016S 
VF "gg V63” 90 
Therefore the 95% confidence limits are 
OR- =ORx exp(—1.96 x ©) 
~ 0.7 x exp(—1.96 x 0.1643) 
hoal, 








BZ 


OR*+ =ORx exp(1.96 x a) 
~ 0.7 x exp(1.96 x 0.1643) 
0 7 
Hence the 95% confidence interval is (0.51, 0.97). 


(b) The observed value of the test statistic for 
McNemar’s test is calculated using Formula (7.5): 
2 2 

(f—a|— 1)? _ (\63-90|-1? _ 4 
fog 63 + 90 
The null distribution of the test statistic is y7(1). 
From ‘Table 3 of the Handbook, the 0.95-quantile of 
x7(1) is 3.84, and the 0.975-quantile is 5.02, so 
0.025 < p < 0.05. 


This provides moderate evidence against the null 
hypothesis of no association. 





(c) There is a negative association between 
heat-related mortality and participation in group 
activities. The odds ratio is 0.70, with 95% confidence 
interval (0.51, 0.97). McNemar’s test provides 
moderate evidence of association (p < 0.05). 





Solution 2.1 


Time plots are introduced in Subsection 1.1 of Book 2, 
and seasonal plots in Subsection 1.2. Stationarity is 
defined in Subsection 11.1. Moving averages are 
discussed in Subsection 4.1. 


(a) The seasonal plot is layered, with lines for later 
years tending to lie above lines for earlier years. There 
is no clear evidence of a common pattern from year to 
year, and hence little evidence of seasonality from the 
plot. This is not surprising, since the data are 
averaged anomalies, that is, averages of temperature 
differences relative to monthly means for 1951-1980. 
Seasonal effects are removed by taking such 
differences. (Note that it is essential to remove 
seasonal effects since data from different climatic zones 
are averaged and, for example, seasonal effects differ in 
the northern and southern hemispheres.) 

















(b) The time plot displays a broadly rising trend, 
with much irregular variation around this trend. 
Because of this trend, the time series is not stationary 
in mean. The size of the irregular fluctuations appears 
to be slightly greater at the beginning of the time 
series than at the end, so the time series may not be 
stationary in variance. 


(c) The higher the order of a moving average, the 
greater is the degree of smoothing produced. Hence 
Figure 2.3(a) is the moving average of order 361, 
Figure 2.3(b) is the moving average of order 61, and 
Figure 2.3(c) is the moving average of order 121. 
Much of the detail has been lost in Figure 2.3(a), 
which appears over-smoothed. In contrast, many small 
fluctuations remain in Figure 2.3(b), which may be 
under-smoothed. Perhaps Figure 2.3(c) represents the 
best compromise between under-smoothing and 
over-smoothing. 

The level of the anomalies fluctuates around a roughly 
constant value between 1880 and 1910, then rises until 
1940 before declining slightly. After 1970, the level 
rises again sharply. 


Solutions to Activities 


Solution 2.2 


The technique of differencing is introduced in 
Subsection 11.2 of Book 2. Significance bounds for the 
correlogram are defined in Subsection 9.2. Selecting an 
ARIMA model is discussed in Subsection 14.2. 


(a) The time series of first differences is stationary in 
mean, so it is not necessary to take second differences. 
Therefore d = 1. The time series of first differences 
appears to be stationary in variance, as the width of 
the fluctuations appears to be constant. There is no 
reason to doubt the assumption of stationarity in 
correlation, so the time series of first differences may 
be assumed to be stationary. 


(b) (i) For a time series with n time points, the 
significance bounds are at +1.96/,/n, so in this case 
they are at 


+1.96/V1277 ~ +0.055. 


(ii) Table 14.1 of Book 2 (which is reproduced in the 
Handbook) is useful when selecting candidate models. 
The correlogram has a single large negative 
autocorrelation at lag 1. At higher lags, the 
autocorrelations are close to zero. The partial 
correlogram appears to tail off to zero. The 
autocorrelation function of an MA(q) model is 
non-zero at lag q and zero at all lags k > q, and the 
PACF tails off to zero in magnitude, so a candidate 
model for the stationary series is MA(1). Thus an 
ARIMA(0, 1,1) model for the series of anomalies is a 
good first choice. 


Solution 2.3 


Fitting and checking an ARIMA model is described in 
Subsection 14.3 of Book 2, and the Ljung—Box test is 
introduced in Subsection 9.2. ARIMA models are 
described in Subsection 14.1. Caution in forecasting is 
advised in Subsection 7.3. 


(a) The sample autocorrelations are close to zero, 
though those at lags 1, 4 and 15 just cross the 
significance bounds. There is also some pattern in the 
correlogram: autocorrelations at lags 2 to 10 are all 
negative. The Ljung—Box test provides strong evidence 
against the null hypothesis that the autocorrelations 
at lags 1 to 20 are all zero. Thus the 1-step ahead 
forecast errors do not constitute white noise. Although 
most of the correlation structure has been accounted 
for by the ARIMA(0, 1,1) model, there is some 
residual autocorrelation that has not been accounted 
for, and hence the model may not be adequate. 


(b) (i) The model formula for the time series of first 
differences Y; is 
Y, — 0.001 = 0.840 (Y;—ı — 0.001) 
+ Z, — 1.458Z4_1 + 0.4682Z;_2, 
where Zą ~ N(0,0.1517). 
(ii) The forecasted anomalies increase steadily to 
December 2050. The prediction limits are wider for 


more distant forecasts, representing increasing 
uncertainty. 


The assumption underlying these forecasts is that the 
processes that give rise to global warming that have 


occurred in the recent past (whatever these processes 
are) will continue unchanged into the future. Thus the 
forecasts do not account for a possible acceleration in 
the processes giving rise to global warming, nor do 
they account for a reduction in such processes. 


Solution 2.4 


Exponential smoothing using Holt’s method is 
described in Subsection 7.1 of Book 2. Obtaining 
predictions for transformed time series is described in 
Subsection 9.3. 


(a) The data represent annual totals, so there is no 
annual seasonality. ‘here is a marked upward trend, 
which is roughly linear. Between 1910 and 1950, the 
level was more variable, perhaps due to the two world 
wars and economic depression. 


(b) Holt’s method is appropriate because there is a 
broadly increasing trend. Simple exponential 
smoothing is inappropriate for the same reason, and 
Holt—Winters is inappropriate because there is no 
seasonal component. 


(c) The starting value for the level is xı = 5.46; the 
starting value for the slope is 


T= pr = 5.49 — 5.46 = 0.03. 





(d) Since a = 1 and y= 0, the level depends entirely 
on the most recent value and the slope depends 
entirely on the starting value. That is, the 1-step 
ahead forecast at time t+ 1 is 


ES =i 0.03. 


(e) The forecasted total emissions of CO2 for 2050 is 
exp(10.3060) ~ 29912 million tons. The key 
assumption is that the underlying increasing trend will 
remain unchanged to 2050. 


Solution 2.5 


Comparing observed and forecasted values of time 
series is discussed in Subsection 14.3 of Book 2. 


(a) The model with natural factors only does not fit 
the data well. In particular, it does not forecast the 
increase in temperature anomalies that has occurred 
since 1970. The model with just anthropogenic factors 
forecasts the increase in temperature since 1970, but 
not the increase that took place in the first half of the 
20th century. The model with both natural and 
anthropogenic factors forecasts both the increase in 
the first half of the 20th century and the increase since 
1970. The forecasts from this model match the 
observations reasonably well. 





(b) Based on current understanding as represented by 
mathematical models, both natural and anthropogenic 
factors appear to contribute to changes in global mean 
temperatures. In particular, anthropogenic factors are 
required to explain the increase in mean temperatures 
since 1970. This suggests that anthropogenic factors 
play a major role in climate change. In fact, other 
work has shown that most of the temperature rise 
associated with anthropogenic factors is due to the 
increase in greenhouse gas emissions. 
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Solution 3.1 


Matrix scatterplots are introduced in Subsection 2.2 of 
Book 3. The correlation matrix is discussed in 
Subsection 4.3. 


(a) There is a positive relationship between winter 
pressures in the Azores and Gibraltar. There are 
negative relationships between winter pressures in 
Iceland and the Azores, and between Iceland and 
Gibraltar. It is unclear whether there are any other 
relationships. The relationships between 
measurements taken in the three locations appear 
stronger in winter than in summer. 


(b) The largest correlation (in absolute value) 
between winter and summer measurements in the 
same location is 0.152, for Iceland. The smallest 
correlation (in absolute value) between winter 
measurements in two locations is 0.642 (between the 
Azores and Gibraltar). In summer, the correlations 
between measurements in two locations range in 
absolute value from 0.170 to 0.357. Thus the 
relationships are stronger between locations (at the 
same time of year) than between seasons (at the same 
location). The relationships are strong in winter and 
much weaker in summer. 


Solution 3.2 


Principal component analysis and standardization are 
discussed in Subsection 7.3 of Book 3. The maximum 
number of principal components that can be found 
and strategies for deciding how many principal 
components to retain are discussed in Subsection 8.3. 
The percentage variance explained by a principal 
component is defined in Subsection 7.2. The 
cumulative percentage variance explained by two 
components is defined in Subsection 8.1, and by more 
than two components in Subsection 8.2. 











(a) (i) For Iceland and the Azores, the mean 
atmospheric pressure is higher in the summer months 
than in the winter months; the converse is true for 
Gibraltar. The mean pressure is lower in Iceland than 
in the other two locations in the winter and in the 
summer. The atmospheric pressure is more variable in 
winter than in summer in all three locations; it varies 
most in Iceland, and least in Gibraltar. 


(ii) The variances of the six variables are very 
different, so it is advisable to base a principal 
component analysis on standardized data. Otherwise, 
the results will be dominated by the Iceland data. 


(iii) The data have dimension p = 6, so at most six 
principal components can be obtained. 


(b) (i) The scree plot displays a rather indistinct 
elbow at component 4. Thus on the basis of the scree 
plot, three components would be retained. 


(ii) When applying Kaiser’s criterion to standardized 
data, only components with variance greater than 1 
are retained. From the scree plot, it is apparent that 
the first three components have variance greater 

than 1, but that component 4 has variance less than 1. 
Therefore three components should be retained. 
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(c) Since the principal component analysis is based 
on standardized data and there are six variables, the 
total variance is 6. The percentage variance explained 
by each component is shown in Table 8.2. 


Table S.2 Percentage variance explained 


Principal Variance PVE (%) 
component 
1 2.368 39.5 
2 1.351 22:0 
3 1.131 18.9 


The cumulative percentage variance explained by the 
three components is given by 


V(¥%1) + VY) + V(Y¥3) 


total variance 


4.85 


~ 80.8%. 


(d) Due to large differences between the variances of 
the six variables, it is advisable to base a principal 
component analysis on standardized data. Three 
principal components explain over 80% of the 
variance, and hence provide a good approximation to 
this data set. 


CPVE = x 100% 





Solution 3.3 


Some guidance on how to interpret principal 
components is given in Subsection 8.2 of Book 3. The 
second principal component is defined in 

Subsection 8.1. 


(a) The largest loadings (in absolute value) for the 
first component are for the three winter 
measurements. The loadings are positive for the 
Azores and Gibraltar, and negative for Iceland. Hence 
this component represents a south minus north 
difference in winter pressures. 


The second component has large loadings (in absolute 
value) for summer pressures in Iceland and the Azores. 
The loading for the Azores is positive, and that for 
Iceland is negative. Therefore this component 
represents a south minus north difference in summer 
pressures. 


The third component has a single large loading, for 
Gibraltar in summer. So this component represents a 
summer Gibraltar effect. 


(b) (i) The first and second principal components 
are uncorrelated by definition, so you would not 
expect to find any linear relationship between them. 


(ii) The point for 1902 corresponds to an average 
value for the first component and a low value for the 
second component. In that year, there was a low 
south—north pressure difference in summer. 


The point for 1990 corresponds to a high value for the 
first component and an average value for the second 
component. Thus in 1990, there was a large 
south—north pressure difference in winter. 


Solutions to Activities 


(c) (i) There is a positive linear relationship between 
the mean January to March temperature in central 
England and the value of the first principal 
component. The point in the top right-hand corner of 
Figure 3.5 corresponds to the year with the largest 
value of the first principal component. From 

Figure 3.4, the year is 1990. 


(ii) The correlation is large and positive, indicating a 
strong positive relationship between the mean January 
to March temperature in central England and the 
value of the first principal component. 


Solution 3.4 


The maximum number of useful discriminant functions 
that can be found is discussed in Subsection 12.2 of 
Book 8. The interpretation of discriminant functions 
based on group-standardized variables is described in 
Subsection 11.3. Calculation of the within-groups 
variance, the between-groups variance and the 
separation is described in Subsection 10.3. 





(a) There are three variables and three groups, so 

p = 3 and G=3. The maximum number of useful 
discriminant functions is the minimum of p and G — 1, 
which is 2 in this case. 


(b) The largest coefficients are for Gibraltar and 
Iceland. Hence this discriminant function represents 
primarily a contrast between pressures in Gibraltar 
and Iceland. 


(c) The grand mean is 
Niti + Negro + N3%3 
N 
46 x (—1.6083) + 46 x 0.0446 + 45 x 1.5984 
7 46 + 46 + 45 


Sl 


~ 0.0000. 


Note that it is not a coincidence that the grand mean 
is 0. This is a consequence of the NAO index being 
defined as the difference between two standardized 
variables. 





The between-groups variance is given by 


G 
1 Sn (T BY 











1 
ecu (46 x (—1.6083)? + 46 x 0.04467 
+ 45 x 1.5984") 
æ 1:7406. 
The within-groups variance is given by 
1 G 
o2o S “(ng — 1)82 
g=1 
l 
= ((46 — 1) x 1.9031 + (46 — 1) x 1.4736 
lof =3 
+ (45 — 1) x 1.4033) 
~ 1.5948. 
Thus the separation achieved by the NAO index is 
Ve 1.7466 
— yx ~ 1.095. 
Vy 1.5948 a 


(d) The total separation is 1.171. Thus the 
percentage separation achieved by D is 


_ 1.170 


PSA l x 100% ~ 99.9%, 





and the percentage separation achieved by the NAO 
index is 


= 1.095 


PSA i x 100% ~ 93.5%. 





Hence the NAO index achieves good separation, but 
not as good as D. It cannot do better than D, since D 
is chosen to maximize the separation. 


Solution 3.5 


The interpretation of within-groups and 
between-groups covariance matrices is discussed in 
Subsection 10.4 of Book 3. One-dimensional 
approximations to bivariate data are discussed in 
Subsections 6.1 and 6.2. 


(a) From Figure 3.7(a), there is a negative linear 
relationship between the annual mean pressures in 
Darwin and in Tahiti. From Figure 3.7(b), there is a 
positive linear relationship between monthly mean 
pressures in these locations. 


Thus the relationship is positive in one case and 
negative in the other. 


(b) The monthly data comprise January means, 
February means, March means, and so on. A 
scatterplot for the January means would show that 
there is a negative relationship between the January 
mean pressures in Darwin and Tahiti. A scatterplot 
for the February means would also show a negative 
relationship. In fact, for each of the twelve months, a 
scatterplot of the monthly means would show a 
negative relationship. However, there is considerable 
variation between months, so when the data for all 
twelve months are plotted together on a single 
diagram (as in Figure 3.7(b)), the within-month 
variation is not evident: the between-month variation 
is dominant, and the relationship is positive. For each 
year, the annual mean is obtained by averaging the 
monthly means for the year. Thus plotting the annual 
means (as in Figure 3.7(a)) reveals the within-month 
variation (averaged over the twelve months). 


(c) The between-month covariance is positive (3.00), 
whereas the within-month covariance is negative 
(—0.32). This confirms the suggestion made in 

part (b) that the within-month and between-month 
relationships between pressures in Darwin and Tahiti 
are in opposite directions. Also, the between-month 
variances are larger than the within-month variances. 
This explains why the between-month variation 
dominates in Figure 3.7(b). 


(d) A good one-dimensional approximation should 
preserve (at least approximately) the ordering from 
the top left-hand corner to the bottom right-hand 
corner of Figure 3.7(a) (or vice versa). This is achieved 
by obtaining the following difference: mean pressure in 
Tahiti minus mean pressure in Darwin. Low values of 
the index correspond to years with low Tahiti 
pressures and high Darwin pressures (that is, the years 
in the top left-hand corner of Figure 3.7(a)), and high 
values of the index correspond to years with high 
Tahiti pressures and low Darwin pressures (the years 
in the bottom right-hand corner of Figure 3.7(a)). 


35 


Statistical perspectives on climate change 


Solution 4.1 


Prior distributions and their interpretation are 
discussed in Subsection 3.1 of Book 4. Normal priors 
are introduced in Subsection 5.1, and the assessment 
of the parameters of a normal prior is described in 
Subsection 6.2. 


(a) The prior for the mean annual rise in sea level 
due to natural factors is less peaked than the prior for 
the mean annual rise due to anthropogenic factors. 
Thus there is more uncertainty about the contribution 
of natural factors than about the contribution of 
anthropogenic factors. For this reason, the stronger 
prior is that relating to anthropogenic factors. 


(b) Experts are rather uncertain about the 
contribution of natural factors to the change in sea 
levels during the period 1910-1990. The most likely 
values for the mean annual rise due to natural factors 
range from —1 mm/yr to +1 mm/yr. The most likely 
value is close to zero (but is negative). 


Experts are much less uncertain about the 
contribution of anthropogenic factors. The most likely 
values lie roughly in the range 0.5 to 1.25 mm/yr, and 
the most likely value is about 0.75 mm/yr. 


(c) (i) First, the overall mean rise in sea levels can, 
in principle, take positive or negative values. Its range 
is therefore unrestricted. Secondly, the 0.25-quantile 
and the 0.75-quantile are located symmetrically about 
the mode. Unrestricted range and symmetry are 
characteristics of the normal distribution. Thus a 
normal prior may be appropriate to represent prior 
beliefs about the mean annual rise in sea levels. 


(ii) The mean a of the assessed prior is equal to the 
mode, so a = 0.70. 


The variance b which matches the assessed quantiles L 
and U is obtained using Formula (6.1) of Book 4: 


2 
p— (= — =) | 
22 

where L is the assessed a/2-quantile and U is the 
assessed (1 — a/2)-quantile of the prior, and z is the 
(1 — a/2)-quantile of N(0,1). In this case, U is the 
assessed 0.75-quantile of the prior, so z is the 
0.75-quantile of N(0,1), which is 0.6745 (see Table 2 
of the Handbook). Since U = 0.99 and L = 0.41, 
Formula (6.1) gives 


2 
0.99 — 0.41 
b= | ———— ] ~ 0.185. 
( 2 x 0.6745 

Thus the assessed prior is N(0.70, 0.185). 





Solution 4.2 


Prior to posterior analyses using the normal/normal 
model are described in Subsection 8.3 of Book 4. 


(a) (i) The posterior for u is given by Formula (8.5) 
of Book 4: 





2 = 2 
udata ~ N (= o-b i 


o2 +nb o2? +nb 
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(ii) Dividing the top and bottom of the expressions 
for the posterior mean and variance by n gives 


osterior mean = (gta + nbz) /n — (m zatiz 
2 2 
posterior variance = _(o%b)/n = (a7 /n) xb 
(o? + nb)/n (o2/n) +b 


(b) Substituting the values % = 0.150, 
o?/n = s*/n = 0.065, and from Activity 4.1, a = 0.70 
and b = 0.185, in the expressions from part (a) gives 
(o7/n) xX a+ bT 

(o? /n) +b 
_ 0.065 x 0.70 + 0.185 x 1.50 
E 0.065 + 0.185 


posterior mean = 


e 1:29. 

(o7/n) xb 
(o? /n) +b 
-0.065 x 0.185 


~ 0.065 + 0.185 
~ 0.048. 


Thus the posterior distribution is N(1.29, 0.048). 


posterior variance = 


Solution 4.3 


Interpreting the results of prior to posterior analyses is 
described in Subsection 9.1 of Book 4. Credible 
intervals are introduced in Subsection 9.2. 





(a) The posterior lies between the prior and the 
likelihood, and is closer to the likelihood than to the 
prior. In particular, the posterior mode (1.29) lies 
between the prior mode (0.70) and the most likely 
value from the data (1.50). 


(b) The posterior is more peaked than the prior, so 
the data have reduced the uncertainty about u. 


(c) The posterior lies well above 0, so the posterior 
probability that u is less than O is very small. Based 
on the posterior, it is extremely unlikely that there 
was a mean annual drop in sea levels. 


(d) (i) The interpretation of the credible interval is 
as follows: with probability 0.95, u lies between 0.86 
and 1.72 mm/yr: 





(ii) Since the posterior distribution is normal, which 
is symmetric and unimodal, the 95% equal-tailed 
credible interval is identical to the HPD interval. 


Solution 4.4 


The interpretation of WinBUGS code is discussed in 
Chapter 8 of Computer Book 4. 


(a) Any quantity in a model definition that is 
followed by the symbol ~ is a stochastic node. 
Therefore the stochastic nodes are x, mun and mua. 
Deterministic nodes and constant nodes are followed 
by the symbol <-. A constant node is a parameter that 
takes a constant value, so the constant nodes are prec, 
precn and preca. There is one deterministic node, mu. 





Solutions to Activities 


(b) In WinBUGS, a normal model is specified in 
terms of its mean and its precision, which is the 
reciprocal of the variance. Hence the priors for u, and 
u, are as follows: 

LU, ~ N(—0.075, 0.162), 


u, ~ N(0.775, 0.023). 





(c) The model may be written as follows: 

X ~ N(y + Mar), 
where X is the sample mean, and the variance 6? 
(which was previously written o7/n) is assumed 


known and equal to 0.065, and the priors for yz, and 
u, are as in the solution to part (b). 


Solution 4.5 


Assessing convergence in MCMC is discussed in 
Section 17 of Book 4. The Monte Carlo standard error 
and the 5% rule of thumb are introduced in Section 11 
and revisited in Section 17. Interpreting estimated 
posterior distributions is discussed in Section 16. 


(a) For each of the parameters u„ and p,, the three 
chains overlap from the start of the simulation. It 
appears that practical convergence is reached very 
quickly after the start of the simulation. Thus a 





burn-in period does not appear strictly necessary. (In 
what follows, no burn-in was used.) 


(b) (i) The 5% rule of thumb states that the MC 
error should be less than 5% of the sample standard 
deviation. In this case, the MC error is less than 1% of 
the estimated standard deviation for both u„ and p,. 
Thus enough simulated values have been obtained. 


(ii) The posterior standard deviations are smaller 
than the prior standard deviations, so the uncertainty 
about u„ and u, has been reduced. ‘The means of u, 
and u, have both increased. The mean of u, has 
changed more than the mean of jz, because the prior 
uncertainty about u„ was greater than that about p,. 


(c) (i) The posterior is roughly symmetric and is 
located almost entirely above zero. The most likely 
values range from 0.5 to 1.2. 


(ii) A 95% credible interval is (0.563, 1.128). This is 
an equal-tailed credible interval. 


(iii) From Table 4.1, the most likely value of the 
mean annual rise in sea levels due to anthropogenic 
factors is 0.85mm. The probability that the mean lies 
in the interval (0.56, 1.13) is 0.95. 
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